rm(list=ls())

# Load packages
library(haven)
library(tidyverse)
library(data.table)
library(lfe)
library(nleqslv)

source("code/config.R")
source(paste0(CODE,"parameters.R"))

df<-read_csv(paste(DATA, "clean_data2.csv", sep=""))

# Solve endogenous variables

## Calculate quantities
temp<- df %>% group_by(sub_area_o) %>% mutate(residential_quant = ((alpha_H_p*R*avg_wage ) + (alpha_H_m*R_w*avg_wage_w ))/residential_rent ) %>% select(sub_area_o,residential_quant) %>% distinct()
df<-left_join(df,temp)

## Extract relevant data
df_d<-df %>% select(sub_area_d,commercial_quant,office_quant,commercial_rent) %>% distinct() %>% group_by(sub_area_d) %>% summarise(commercial_sum=sum(commercial_quant*commercial_rent),office_sum=sum(office_quant*commercial_rent),commercial_quant=sum(commercial_quant),office_quant=sum(office_quant))

## Calculate N, T variables
df_d$share_N_p<- ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_p)/( (1-beta_N)*beta_T*beta_T_p) )/(1+ ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_p)/( (1-beta_N)*beta_T*beta_T_p) ))
df_d$share_T_p<-1-df_d$share_N_p
df_d$share_N_m<-((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_m)/( (1-beta_N)*beta_T*beta_T_m) )/(1+ ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_m)/( (1-beta_N)*beta_T*beta_T_m) ))
df_d$share_T_m<-1-df_d$share_N_m
df_d<- df_d %>% select(sub_area_d,share_N_p,share_T_p,share_N_m,share_T_m)

## Merge data
df<-left_join(df,df_d)

## Save
write_csv(df,paste(DATA, "clean_data3.csv", sep=""))
